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CHAPTER 1 
INTRODUCTION 

In many hydraulic engineering applications, it is required 
to design open channels for super c r i t i cal flow 
condi t ions . Spi 1 lway chutes sewer systems and outlet works are 
some common examples where the flow is 

supercritical .supercritical flow is also encountered in road 
drainage works. Some times an increase in the capacity of hydel 
power plants requires higher velocities of water in the 
conveying structure and this makes flow super criti cal .These 
high velocity flows can also occur in natural channels. The 
flow in steep mountainous streams and in rivers could be 
supercritical during periods of high flows. A typical example 
of this is dam break flow. 

The analysis of supercritical flow is more complicated than 
the analysis of subcritical flow because of a number of 
reasons. Roll waves or slug flow may be present due to 
instabilities which develop when the Froude number is greater 
than a critical value. Supe r cr i t i cal flows are also be 
accompanied by air entrainment and cavitation. The most 
importantcharacteristi cof supercr iti calf lows is the presenceof 
standing waves and large surface disturbances. Many times it 
is necessary to provide lateral change of the boundary 
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confining the flow at supercritical velocities. These local 
changes either in cross section, slope or wall direction of 
the channel results in the formation of standing waves which 
will cross and recross the channel affecting the down stream 
flow significantly. Supe r c r it i cal flow through channel 
junctions is a good example ofthe above phenomenon, 
cross-waves similar to those as in contractions develop and 
pertain a long distance downstream of the channel. Depending 
on the junction geometry and inflow Froude numbers a 
hydraulic jump may form in one of the branches. The side walls 
have to be considerably higher than those normally provided 
because the maximum height of the free surface may be 10 to 20 
times as large as the in flowing depths. Difficulties are 
encountered at junctions in tunnel spillways or closed 
channels due to pulsations, air entrainment, transitions to 
pressurized flow and moving hydraulic jumps, if the height of 
the downstream channel is insufficient. The present study is 
an attempt to develop a mathematical model for the analysis of 
supercr iti cal flows in channel junctions. 

1.1 Experimental and Analytical studies on supercritical 
flows 

The first contribution on the subject of super criti cal flow 
in non prismatic channels originated from investigations in 
the field of gas dynamics. 7Raibouchinsky and Prandtl were the 
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first ones to point out the analogy between supersonic flow of 
air and supercr iti cal flow of water for two dimensional motion 
(Ippen 1953). Prandtl (1953) developed a graphical method for 
the analysis of supersonic flow in gases. Preiswerk (1940) 
extended this method which is based on the character isti cs 
theory to supercr iti cal flow problems in open channels. 
Laitone (1953) critically examinedtheconceptof hydraulic 
analogy and illustrated some of its limitations. Loh (1959)has 
correctly pointed out that the hydraulic analogy is based on 
the assumption of flows withnegl igiblevertical accelerations, 
and this is not valid when the shocks are strong. The 

limitation of the analogy arises f romthef actthatf or 

gases, Euler equations are moreor lessexactmathematical 

representation of the phys i cal phenomenon , but for open 
channel flows, the shallow water equations, in certain 

ci r cumstances , a some what crude approximation. However, the 
experiments conducted by Ippen and Harleman (1950, 1956) 

indicated that the hydraulic theory based on gas dynamics 

analogy works well for a wide range of Froude numbers. They 

concluded that the limitation comes from the fact that no 

information can be obtained about short waves, rollers and 
development lengths using the gas dynamics analogy. Engelund 
and Munch-Petersen (1953) used potential flow theory to show 
that the expression for the angular position of a small 
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stationary wave as derived using the ^as dynamics analogy is 
true for an infinitely wide channel. It is actually afunction 
of Froude number and the width to depth ratio. 

An excellent summary of the early work on supe r c r i t i cal 
flows is presented in the"proceedings of a symposium on high 
velocity flow" (Ippen et al . 1951) . Ippen andKnapp( 1936, 1938) 
and Knapp(1951> studied the super c r i t i cal flow in the bends 
Ippen and Dawson (1951) , Harrison (1966) and Sturm (1985) 
gave design considerations for supercritical flow in 
contraction. Rouse (1938) and Rouse et al (1951) analysed 

super critical expansions and developed a special geometry for 

* 

expansions in order to prevent flow separation and to obtain a 
water profile with minimum of surface disturbances . 

Most of the available work on supercritical flows is 
concerned with the analysis for channel transitions andbends. 
Very little work has been reported regarding the supercr iti cal 
flow characteristi cs in channel junctions . Bowers (1950) 
seems to be the first to have analyzed supercritical flow 
through channel junctions . He found that the formation of a 
hydraulic jump in one of the branches depends on the junction 
geometry and the inflow froude numbers .Schnitter et al (1955) 
conducted experiments in channel junctions where the angle of 
confluence is almost zero, but the lateral channel had a steep 
slope. They observed considerable cross wave pattern in the 
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downstream branch. This wave became strong when the twc 
upstream branches have different in flowing Froude numbers. 
They also gave guide lines for reducing theperturbat ions . 
Behlke and Pritchett ( 1966) analysed supercritical flow ir 
junctions for different angles of confluence and for a range 
of Froude numbers. They presented a theoretical approach for 
determining wave angles in simple junctions based or 
relationships established by Ippen and Knapp(1936). This lee 
to identification of wall pile up zones. Greated (1968) alst 
presented similar relationships for the main wave angle at i 
simple junction. Gerodetti (1978) conducted experiments or 
supercritical flows as encountered in road drainage channels, 
He used two semi-ci r cular channels of top width 15 cm, bottor 
slope of 28%, and an angle of confluence <5=37°. Recently 
Hager (1989) used both theoretical and experimental approaches 
to analyze the main flow features in junctions witl 
thoroughly super cr it i cal flows. The junction is character izei 
by rectangular channels of equal branch width and sharp edge' 
corners. Junctions angles of 22.5°and 45°are cons idered .Thei 
analytical approach , although elegent unfortunately has a 
error. The correct derivation of equations for estimatin 
maximum height of water level in channel junction based o 
Hager’s approach is presented later. Problem of interest i 
open channel flows is the division of flow into branc 



1.2 Analytical approach for naximua, watar l.val in 
channel Junctions. 

The analytical model given below closely follows the 
derivation presented by Hager M989) .Unfortunately has an error 
This is corrected in the present analysis. Derivation is for 


"simple junction" that 

is r a junction 

where 

the two 

inflow 

branches have a common 

invert 

e 1 e vat ion 

and 

there 

are no 

baffles lor piers to deflect the 

f 1 ows 

at 

the 

j un ct i on 

. Also 

the main channel is 

straight 

and 

i s 

of 

constant 

width. 

Streamline curvature effects are 

not 

incorporated. 
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Refering to Fig. 1.1, Upstream of point A, two plane flows 

are present which are characterised by h , F and h,,F,h=flow 
f u u 11 

depth, F=froude number and subscripts u and 1 refer to the 
values in the main and lateral channels upstream of the 
junction respectively. At point A, the two undisturbed flows 
intercept each other, thereby creating a main wave front at an 
angle 9 relative to the main direction. The two incoming 
flows are thus deflected by the angle 9 (U/s branch) and (<5 - 

9) (lateral branch) thereby forming two new wave fronts. Due 
to main wave front at an angle 9 the upstream flow turns once 
the wave front at an angle ft*, is passed and this latter wave 
front hits the walls. It is reflection leaves at an angle fte 
and streamline become parallel to the wall once this second 
wave front is passed. Thus three different domain of flow may 
be distinguished as follows 

(1) upstream zone (index "u" ) 

(2) Zone between fronts with angles 9 and ft*. <index "1") 

(3) Zone between wall and wave front with the angle 



proposed by Ippen and Knapp. (1936) This is illustrated below. 



SUPERCRITICAL, FLOW PAST A V/ALL DE FLECTION 


Fig. 1.2. show the stream line pattern for supercritical 
flow past a wall deflection. Here “the angle by which the 
wall is deflected i.e. the angle by which the flow changes it 
directions and is the resultant wave angle. Subscripts ± and 
2 refer to the values upstream and downstream of the wave 
angle front respectively. The following relationship can be 
used to determine h2 , F 2 and ft if hi,Ri., and w are given. 
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Fig. 1.2. show the stream line pattern for supercritical 
flow past a wall deflection. Here w “=the angle by which the 
wall is deflected i.e. the angle by which the flow changes it 
directions and is the resultant wave angle. Subscripts i and 
2 refer to the values upstream and downstream of the wave 
angle front respectively. The following relationship can be 
used to determine ha , F 2 and ft if hi ,R , and u are given. 


ha/hi = 1/2*C<1 + 8*F* a *Sin 2 /3 )*'' 2 - 13 (1.1) 

h 2 /ht = tary9 / tan(/3 -ta ) 1 .2) 

Fa 2 = y"** CF 4 2 - ( 1/2*y >*<y~1 >*<y+1 )* 3 <1.3) 


These three equations constitute the governing system for 
the determination of the sought quantities ha, Fz , ft for given 
ht,Fi and < 0 . As can be noticed, determination of hi ft and ft is 
implicit Hager and Bretz (1987) derived the following explicit 
approximate equations from the Eqs. 1. 1-1.3 for the conditions 


that ft <45° and Fi >2. 

ha /hi = 42 * Fi*siry? - 1/2 <1.4) 

ft = oi + 3/(242 * Fi) (1.3) 

Fa® = Fi*Cos 2 f3/<4 2*sirv5*< 1 -2.0*/* >) (1.6) 

where 

ft. = 4 2*R.*siryJ ...........(1.7) 


Applying eqs. 1.4-1. 7 to regions "u" and "1" in fig. 1.1, We 
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get 

fa + 3/(2<|2*Fu) ..(1.8) 

yst = 4 SsFufcSirf?! -0.5 (1.9) 

and 

Fa* = Fu*Cos 2 /?* /(«| 2*Sir*3t *C1 - (Sl^t^Ri)" 1 ) (1.10) 

A similar application for regions 1 and 2 gives 

fa = e + 3/(2<| 2*Ft) (1.11) 

yz = -|2#Fi#Sinf3a -0.5 (1.12 

and 

Fa 2 = Ft #C o s Z /?2 / (-J 2#S i r\ftz -s-C 1 - ( 2j 2fia. *F* f 4 )...........( 1 . 13) 

Eqs. 1.8-1.13 can be used to determine h2 and F 2 Once the 
value of &, the angle by which the main flow changes it 
direction due to interaction with the lateral flow is known. 

& depends on the inflow parameters ^'u'^l'^u'^] anc * ^he 
confluence angle 6 .Greated( 1968) derived the following 
equations for simple channel junctions based on the 
assumptions of uniform velocity and hydrostatic pressure 
distribution. 

1~Y 2 » 2#Y 2 *Fu 2 *Sin 2 <9 - 2*f l 2 *Sin 2 *(<5 - S) (1.14) 

where 

Y « hu/hl (1.15) 

Hager (1989) further simplified eq. 1.14 and derived the 

following. 

e/6 = C 1 + (Vu/Vl )#(hu/hl) 1/S D" 1 (1.16) 
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Equation 1.16 is valid provided 6 £ 45° and 
1/(11 - Y 2 ! > iy,z *Fl > 1/S 

In this study Eqs . 1.16 ,1.8 -1.12 are used for determining h 2 
for given inflow conditions and the confluence angle. 

1.4 Numerical Modelling of Supercritical flows 

In spite of the degree of refinement reached in the 
computation of open channel flows and the significant advances 
on the subject, it appears that in many aspects the 

supercr iti cal flows still defy reliable analysis (Lai 1986). 
Steady two-dimensional shallow water equations constitute a 
set of hyperbolic partial differential equations. Therefore, 
character isti c method can be used for their solution. Bagge 
and Herbich (1967), Herbich and Walsh (1972s), Villagas (1976) 
and Daksh inamoor thy (1973, 1977) used this method for 

analyzing steady supercritical flows. Except Dakshinamoorthy 
(1973), all the above investigators were primarily interested 
in super criti cal flows in channel transitions. Dakshinamoorthy 
(1973) makes a reference to the problem of flow division in a 
junction when the upstream flow is supercritical. However, he 
did not analyze the shock wave pattern in junctions where 
super critical flow from one channel interacts with the 
super critical flow from another channel. Character istic method 
requires many interpolations which could seriously affect the 
accuracy of the solution. Also , characteristic methods are 
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not able to compute oblique jumps properly. 

The analogy between gas dynamics and shallow water flow has 
prompted many researchers to borrow techniques developed in 
the former area and apply them to supercritical flows in open 
channels. There are two types of techniques to compute flow 
fields with shocks! "shock-fitting" and "shock-capturing" 
techniques Shock-fitting methods treat all shocks as 
discontinuities and explicitly compute their motion and 
strength. In shock-capturing techniques, shocks are predicted 
as part of the solution and need no special treatment .Ligget 
and Vasudev (1965) are the first to numerically integrate the 
steady two-dimensional shallow water flow equations. They used 
the methods developed by Stoker (1957) and analyzed 
supercritical flow in expansions including the effects of 
bottom slope and friction. Their method, like the others 
mentioned above, can not compute flows with shocks or standing 
jumps. Pandolf i ( 1975) applied shock-fitting methods for the 
computation of steady two-dimensional flow around a blunted 
obstacle in a supercritical flow stream. Demuren (1982) used 
Patankar and Spalding’s method for simulating both 
subcritical and super cr i ti cal flows in open channels. However, 
the ability of numerical scheme to handle discontinuities is 
not clearly demonstrated. Recently, Jimenez (1987) and Jimenez 
and Chaudhry (1988) used a shock capturing technique to 
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simulate supercritical flow in channel contractions and 
expansions. Their technique is able to simulate the shock 
proper ly. Steady supercriti cal flow equations are hyperbolic in 
nature and shock-capturing, marching techniques may be used to 
simulate the oblique jumps or shocks. It can not be used if 
the flow is subcritical in part of the channel since the 
governing equations become elliptic. Therefore, all the 
methods mentioned earlier are applicable only when the flow is 
supercr iti cal in the entire length of the channel. In such 
cases, an unsteady flow model may be used to obtain steady 
flow solutions by using the false transient method. Such 
method is capable of simulating mixed sub and supercritical 
flows if a proper numerical technique is used. In the false 
transient method, some arbitrary initial conditions are 
assumed and the boundary conditions are set equal to the 
initial steady state conditions. Using the governing equations 
for unsteady flow, the system conditions are computed for 
sufficient length of time until the variation of flow 
conditions is negligible and the conditions have converged to 
steady values corresponding to the boundary ponditions. Many 
numerical schemes are available for solving two dimensional 
unsteady flows( Katopodes 1984, Garcia and Kahawita 1986, 
Fennema and Chaudhry 1990). Recently, Bhallamudi and Chaudhry 
(1991) -<sed this approach -for computing steady supercritical 
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flows in channel expansions and contractions » Their results 
matched satisfactorily with the results obtained using a 
steady flow model. The ability of the model to simulate mixed 
sub and supercritical flows was demonstrated by simulating 
hydraulic jump in an expansion. 

From the above discussion it can be summarized that most of 
the earlier experimental and analytical work in supercritical 
flows has dealt with channel transitions and bends. There are 
some recent works on super c r i t i cal flow junctions (Hager 
1989). However, analytical methods are based on many 
assumptions and do not describe the complete two dimensional 
flow characteristi cs in channel junctions. Therefore, 
numerical methods should be used for the analysis. Almost all 
of the earlier numerical solutions are for channel 
transitions. To the authors* knowledge, no attempt has been 
made as of yet to numerically simulate supercr iti cal flow in 
channel junctions. 

1.5 Present Study 

The aim of this study is to develop a mathematical model 
for the computation of steady two dimensional supercritical 
flows in channel junctions.lt is assumed that the shallow 
water flow equations describe these flows, although it is 
recognized that this assumption may restrict the use of 
mathematical model. Therefore, one objective of this study is 





CHAPTER 2 
GOVERNING EQUATION 

The shallow water equations describe the gradually varied 
flow in open channels. They follow from the application of 
the principles of the conservation of mass and momentum. 
Since supercritical flow in a junction in open channels is 
characterized by relatively large changes in depth over short 
distance, shallow water equations are not exactly applicable. 
However, solution of more general equations describing flows 
in open channels like Boussinesq type equations which include 
vertical acceleration ef f e ct (neg le cted in gradually varied 
flow equation) is not simple even for one-dimensional flows 
(Gharangik 1988, Abbot et al. 1978). It becomes a lot more 
complicated for two-dimensional case. Therefore, Shallow 
water equations are used in this study as a first step 
towards analyzing supercritical flows in junctions. The limit 
of applicability is established in an indirect way by 
comparing numerical results with the experimental data. 

In this chapter, the shallow water equations and the 
assumptions on which they are based are presented. This is 
followed by the discussion on the limitations of the shallow 
water theory. The coordinate transformation procedure for 
transforming the physical plane into a rectangular 
computational plane is also presented. 
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-gh< S - S. ) 

= . * ox f x 

-g h (S 0y - S fy ) 

Subscripts t, x and y in eg .2. 1 . indi cate the partial 
derivatives. In the above equation; 
t= time 

u r v = depth averaged flow velocities in the x,y directions 
respectively. 

h = water depth measured vertically, 
g = acceleration due to gravity. 

^ox'^oy = Channel bottom slopesalong longitudinal and 
lateral directions respectively. 


Channel slopes S and S are given by 
h ox oy 3 


S = Sinot (2.2) 

ox x 

S = S i not (2.3) 

oy y 

In which 

( x , y ) »« y ( x , y ) = angle between the bottom of the channel 
and the (x,y> direction along longitudinal and 
lateral direction respectively. (See Fig. 2.1) 


Friction slopes are estimated using any empirical equation 



like Chezy* s eqution. 


S f = <uj u”+v a )/c a h (2.4) 

S r = (vj U 2 +V 2 )/c a h (2.5) 

fy 

Where 

c = Chezy's coefficient. 

S^ i<r S^y= friction slopes in the longitudinal and lateral 
directions respectively calculated from chezy type equations. 

2. 2 Assumptions 

Equation 2.1 is based on the following main assumptions 

1) The accelerations in the vertical direction are 
negligible, making the pressure distribution in vertical 
direction hydrostatic. 

2) . Velocity distribution is uniform along the flow depth. 

3) The fluid is incompressible and has homogeneous density. 

4) The channel bottom slope is rigid and small. 

5) Only shear stresses due to horizontal velocity 
components are significant, other viscous effects are 
neglected. 

6 ) The frictional resistance is approximated by the chezy 
or the Manning equations. 

A discussion on the fore mentioned assumptions is being 



presented in the following paragraphs : 


2.2.1 Hydrostatic pr«s a ur« distribution! Deviation from 
hydrostatic pressure distribution is negligible as long as 
regions of high curvature are not present. Ligget (1975) 
showed that the above assumption is valid as long as a 
"shallowness" parameter ha/lo(ho = water depth and lo =a 

characteristic length.) is small. However, problem arises in 
defining the characteristic value in each case. Jimenez( 1987) 
estimated the errors involved in supercritical flow 

computations if shallow water theory is used. He compared the 
wave angles obtained using shallow water theory to those 
obtained using the more regorous theory of Engelund and 

Munch-Pete rsen ( 1953 ) based on Potential flow assumption'. It 
was concluded that the shallow water theory reasonably 

represents the supercritical flow if depth to width ratio is 
of the order of 0.1 and the Froude number is not close to one. 
The error is of the order of 20'/. if the depth to width ratio 

is increased to 0.2. This error is manifested in the wave 

length of the resulting wave patterns. Ippen and Harleman 
(1956) also showed that this theory provides accurate results 
except in the vicinity of the slope itself . Some details are 
lost, but the overall results are adequate for engineering 
purposes (Cunge 1975). One of the main aims in the design of 
hydraulic structure having super c r i t i cal flow is to reduce 
shocks in water surface profile as much as possible. The 



methods based on shallow water equations provide useful 
results for this purpose. 

2.2.2 Depth averaging and -turbulent effecti The depth 
integration of the three-dimensional equations of motion 
produce higher order terms known as effective stresses which 
are not included in equation 2.1. In a detailed examination of 
this problem, Flokstra (1977) identified three types of 
"effective stresses" that arise from the averaging process. 

1) Laminar viscous stresses 

2) Turbulent stresses 

3) Stresses due to depth averaging of 
advective terms. While the first type of effective stresse can 
be neglected, the other two are important. Consideration of 
turbulent stresses requires the use of turbulence closure 
models which express stresses as functions of the mean flow 
variables. Rastogi and Rodi (1978), Puri and Kao (1984) and 
Leschziner and Rodi (1979) used "K-*>" models for describing 
turbulence. Vreugdenhil and Wijbenga (1982) used a much 
simpler constant eddy viscosity concept. However, the above 
closure models are tested only for subcritical flows and no 
information is available for supercritical flows. 
Non-inclusion of turbulent terms leads to omission of side 
wall boundary layer effects. Ippen and Harleman (1930) 
analyzed the case of oblique jump and concluded that reduction 



of velocity across the boundary layer reduces the vertical 
ac ce 1 e rat i ons at the point of wall deflection. This makes the 
shock front steeper than in the case when the boundary layer 
is not included. 

Ippen and Harleman (1950) studied the effect of non-uniform 
velocity distribution also. The measurements of velocity 
across the depth on the front and back side of the jump 
indicated that the momentum correction factor ranged from 
1.007 at F = 6.3 to 1.015 at F =2.0. Therefore, they concluded 
that the effect of velocity distribution may be considered 
negligible. However, it should be noted that although the 
effect of magnitude variation of velocity across the depth is 
not significant, the variation in the direction of resultant 
velocity may have some effect. This is especially true for 
flows in channel junctions as found in the experiments by 
Hager (1989). Fig 4.13 shows the velocity vectors at different 
depths for super cr iti cal flow in a channel junction. As can be 
seen, the velocity is almost parallel to the side wall near 
the channel bottom and it is pointing towards the wall near 
the surface. Equation 2.1 can not simulate this effect and the 
numerical results will be in error. As mentioned earlier, one 
of the purposes of this study is to partially estimate the 
errors due to above assumptions by comparing numerical results 
with the experimental data. 



Before this discussion is closed, it should be pointed out 
it is not possible to simulate separation if the 
.effective stresses are not included (Flokstra 1977) Any 
circulating pattern in numerical results obtained without the 
inclusion of the effective stress terms is only an indication 
□f numerical errors. It should also be mentioned that only 
channels with vertical side walls are considered in the 
present study. There are several complications in the case of 
non rectangular channels from a computational point of view 
because of the existence of zones that become dry or wet and 
theoretically because the depth may become zero at the edge of 
the channel. 

3.3 Coordinate transformation 

In this study , a finite-difference method is used to 
solve eq. H.l. Physical boundary as well as the finite 

difference grid in cartesian coordinate system for a channel 
junction are shown in Fig.H.3. 



We notice that it becomes necessary to represent the 
physical boundary in angular portion in a stair step fashion . 
It leaves sharp corners (marked 1 on Right side and 2 on Left 
side). Sharp corners introduce numerical disturbances; which 
spread quickly to entire mesh. This may result in numerical 
instability (Roache 1972). This is particularly so for 
super cr i t i cal flows which are hyperbolic in nature . An 
alternative method of arranging the grid could be as in 
Fig. 2. 4 ' 



This a|gain leaves wedges at the junctions of two channels 
and also requires different coordinate directions for straight 
and angul|ar channels. Aforementioned problems can be avoided 
by using [a coordinate system such that the coordinate axes 
coincide with the boundar i es ( Ancle r son et. a l . 1984). Use of a 
s imp 1 e algebra! c coordinate transfo r mat i on suf f i ee s i n the 
present case. 

The Tol lowing transformation of independent variables 
onverts the original physical plane ( E . 5a) into a rectangular 
ornputat, i ona 1 plane (Fig E.5 b. ) 






The governing equations must be expressed in terms of the new 
independent variables for their application. This is done by 
applying the chain rule of partial dif f erentiation as follows! 


d/9« = <o * a /at; ) + (rpt #s/&n) (£.8) 

d/dy = <Cy *e/dt; ) + (r) y * d /&r) ) (2.9) 


Subscripts denote partial de r i vat i ves -Eqs . 2.6 and 2.7 give 


C* = 1 » = 0 (2.10) 

Cy = Cota f r)y = 1/Sina ...(2.11) 

a/4x = S/at: ;^/dy = Cota*0/««; + ( 1/Sinot )*d/&n (2.12) 


Application of equation 2.12 to equation 2.1 and 

subsequent rearrangement of terms yields the following 
conservation equation. 

E + F + G + S = 0 (2 . 13 ) 

1 x y 

In which 



„ uh + hv*Cota 
u n + gh*/2 +huv#Cota 

n O 

uvh + (hv + gh /2)Co tot 


G 


vh/Sinct 

uvh/Sin« 

( v 2 h + gh 2 /a)/Sina 



27 , 


S 


-gh (S - S. ) 
a ox f :< I 

-gh (S - S. ) 
d oy fy -* 


Equation 2.13 is solved in this 
finite-difference method. The details of the 
are presented in the following chapter. 


study using a. 
numerical scheme 



CHAPTER 3 


NUMERICAL SCHEME 

The shallow water flow equations presented in the previous 
chapter constitute a set of hyperbolic partial differential 
equations. Analytical solutions for these are available only 
for highly simplified cases. Therefore, they are usually 
solved by suitable numerical techniques. Either 
finite-difference or finite-element techniques can be used for 
this purpose. Finite element techniques for open channel flows 
are only in an infancy stage and their application is 
complicated, especially for flows with shocks such as those 
considered here. Therefore, a finite-difference method is 
chosen in this study. We are interested only in the steady 
supercritical flow characteristi cs in junctions. These steady 
flow characteristics are obtained by the false transient 
method. Any initial conditions, such as constant velocity and 
a constant depth, through out the system are assumed and the 
boundary conditions are set equal to the steady-state 
conditions. Using the two-dimensional unsteady flow equations, 
the system conditions are computed for sufficient length of 
time until the variation of flow conditions is negligible and 
the conditions have converged to the steady values 
corresponding to the boundary conditions. In other words, the 
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time is used as an iteration parameter. 

The supercritical flow in channel junctions is associated 
with a complex pattern of shocks. Therefore, it becomes 
convenient to use a shock capturing technique rather than a 
character ist i c method with shock fitting procedures. According 
to the theory upon which these methods are based (Lax 1954, 
Lax and Wendroff 1960, Abbott 1974) three steps should be 
considered. 


1) The governing partial differential equations 


should be in conservation form. 

Reason for 

this 

is 

that 

Eq.2.1. expresses the 

laws of 

conservation 

of 

mass 

and 

momentum, and these are 

the quantities which 

are 

actually 


conserved through a jump. 

2) . The numerical scheme should be conservative. In 
other words, the scheme should exactly maintain the 
discretized version of conservation laws for any mesh size and 
over any arbitrary finite regions. 

3) The numerical scheme should be dissipative in 
nature in order to represent the energy losses associated with 
hydraulic shocks (Roache 1972). This dissipation should be 
very small in regions of gradual variation, but comparati vely 
large near a discontinuity (Abbot 1975). Some numerical 
schemes, particularly the first order schemes, have "inbuilt" 
dissipation mechanism. The schemes which do not have 



inbuilt" dissipation require explicit inclusion of viscosity, 
known as "artificial viscosity" for successful application. 

Shock capturing techniques can not precisely resolve the 
shock location and they have a tendency to produce spurious 
oscillations near the strong shocks. The first objection is 
valid when the shocks are true discontinuities. However, this 
is not the case in open-channel applications. The spurious 
oscillation near a shock front originate because the 
conservative variables being differentiated across the shock may 
have discontinues de r i vat i ves ( Kut 1 e r 1975). Generally, these 
osci l\ations are a characteristic of second and higher order 
schemes.: These oscillations can be reduced by 
1) introducing artificial viscosity and 

S) aligning one of the coordinate axes in the direction 
of shock ( Anderson et. al . 1984). 

Application of the finite-difference methods require the 
use of a finite difference grid as shown in Fig. 3.1. 
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The values of the dependent variables at a grid point are 
identified by the node (i,j). i indicates the position along 
< direction and j indicates the position along y) direction. AC 
and An are mesh step sizes in the C and Y) directions 
respectively. The dependent variable are known at any time 
either as initial conditions or through computations for a 
previous time step. The values of the variables after a time 
step of At i.e. at time t + At are determined using the finite 
difference technique as described in the next section. 

3. 1 Mac Cor mack Scheme 

A large number of numerical schemes are available to solve 
hyperbolic system of equations. In general, the choice of a 
particular method is directed by several factors such as the 
simplicity of programing, accuracy, efficiency and the 
particular application under consideration. Besides these 
considerations accuracy, convergence, diffusive properties, 
stability etc. of a numerical scheme are also important 
consideration. MacCommark scheme is used for the problem under 
study. It is an explicit, second order accurate shock 
capturing scheme for solving hyperbolic equations. This scheme 
is explained in detail elsewhere and only a brief description 
is given here. It comprises of a predictor and a corrector 
step. In the predictor part, backward finite-differences are 
used for approximating spatial partial differential terms. 
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In the corrector part, they are approximated by the forward 
difference scheme, using predicted values. With reference to 
the finite-difference grid shown in Fig. 3.1 these 
approximations are 


Predictor part 

*//*£ = </ k (i, j> - / k a-4,j ) ) / AC (3.1) 

df/dv) = ~ / k a,j-i ) ) / Af) (3.2) 

Sf/dt = (/ (i,j ) - f k <i,j))/ At (3.3) 


where / is used for any function of which the partial 
derivative is required. The subscript k denotes the value of 
the function at known time level, i.e. at time t , subscript p 
refers to the values at the end of the predictor step. At is 
the computational time step. 

Corrector part 


*//*C = (/ (i+M ) -/ <t„j))/A£ (3.4) 

P P 

. . tf/&n = </ (m** ) -y <L,j)>/A n (3.5) 

p p 

Sf/St = (/ (t,j> - / (i.j))/ At (3.6) 

G P 


Where subscript e indicates the value of the function at the 


end of 

corrector step. The 

value of 

the function 

at 

the 

unknown 

t ime level i.e 

at time 

t+At is given 

by 

the 


following 
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/ k+i = 0.5*(/ k +/ c ) (3.7) 

Substitution of Eqs.3.1 to 3.7 in Eq. 2.1 leads to explicit 
algebraic expressions for the determination of h, u and v at 
the unknown time level t + At if the values are known at time 
t. This requires that initial conditions be given for 
starting the unsteady computations. Once the values are known 
at t+ At these are used to determine the values at time 
t+2#At. This marching in time is continued till steady values 
of h, u and v are obtained. 

The above procedure can be used for all nodes i=2, max(i) 
and j=H, max ( j ) where max (i) is the last but one node in £ 
direction and max (j) is the last but node in r) direction. 
Boundary Conditions have to applied at the boundary nodes. 
Initial and boundary conditions are explained in the next 
section. 

3. 2 Initial and Boundary Conditions 

As mentioned earlier, initial conditions have to be 
specified in order to start the computation. Since we are 
using a false transient method to obtain steady state flow 
characteristics, any initial conditions such as constant 
velocity and constant depth equal to the boundary value may be 
assumed for all nodes. 

The determination of u, v and h at the boundary nodes 
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during unsteady computation involves the application of 
boundary conditions. The treatment of the boundaries is one of 
the most important aspects in the successful application of 
any numerical technique (Moretti 1969). Hyperbolic equations 
are particularly sensitive because errors introduced at the 
boundaries are propagated and reflected throughout the grid, 
and in many cases instabilities may result (Anderson et 
al.1984) The literature covering the numerical applications to 
hydraulic problems refers to this subject rather scarcely. Two 
types of boundaries are encountered in the present study. They 
are 

1) Flow boundaries and 

2) Solid side wall boundaries. 

The treatment of the boundaries in the numerical mode is 
described below. 

з. 2. 1 Flow Boundaries: 

Stoker (1957), Abbott (1979) and Cunge et. al. 
(1980) discussed the mathematical basis for the specification 
of open boundaries in shallow water wave applications. Based 
on the characteristic theory, it can be shown that for two 
dimensional supercritical flows, three boundary conditions 
have to be specified at the inflow boundary and none at the 
out flow boundary. 

и, v and h are specified at the inflow sections of main and 
lateral channels as upstream boundary conditions, u, v and h 
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at the downstream end are obtained simply by extrapolation 
from the values at the interior nodes already calculated by 
MacCormack scheme. Although, this procedure is not exact, the 
errors introduced are limited to only few nodes at the 
downstream end. This is because the flow is super critical and 
has an upstream control. A disturbance can travel only in the 
downstream direction and not in the upstream direction. 

3.2.2 Solid side wall boundary 

The basic requirement for a solid boundary is that 
there is no flow normal to it. This is expressed as 

tanS = v/u .....(3.8) 

where B is the angle between the wall and the axis. This is 
called the slip condition and is the proper boundary condition 
at a side wall since all the stresses other than the bottom 
stress are neglected. For inviscid gas flow applications, 
several wall boundary techniques have been developed Abbett 
1971, Moretti 1969, Roache 1972, and Anderson et al . 1984). 

All these techniques try to enforce in one way or another the 
basic requirements given by eq 3.8. The problem 

arises because, in order to apply the numerical 


scheme 

at 

the 

grid 

points 

along 

the wall, the 

values 

of 

all 

the 

variables 

are 

required and the 



Therefore, it 


equation does not provide enough information. Therefore, it is 
necessary to find a way to compute these values using 
information from the interior points plus the ' boundary 
condition. The problem was studied by Abbett (1971) who compared 
and classified them into four broad categories as follows; 

1) Reflection or image point procedures. 

H) One sided derivatives procedures , 

3) Procedures based on the method of characteristics and 

4) Abbet's procedure. 



Jimenez (1987) tested several techniques for the case of 
steady supercritical flow in transitions and concluded that the 
Abbe t * s procedure is best suited for this purpose- However y this 


! 



procedure is • appl i cable only For steady Flow conditions and is 

not suitable for unsteady Flow computation done here. The 

reflect ion procedure seems to be more suitable than the other 

techniques and is used herein. In this procedure, all 

non-conservative flow variables other than the normal velocity 

component (i e u, and h) are considered even functions with 

respect to the wall. The normal velocity is assumed an odd 

function to ensure zero normal velocity at the wall. In 

implementing reflection, an imaginary row of points is 

considered as shown in fig. 3.2. values of the variable at the 

$ 

imaginary points are «' • ■ -I according to the above symmetry 

conditions. For the purpose of illustration let us consider 
boundary next to J=1 row (Fig. 3. 3.) 



Since wall under consideration is aligned with the X axis, 
boundary condition is given by v = O at the boundary . The 
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reflection procedure is slightly modified here because of 
transformations of independent coordinates. The mathematical 
representation of reflection procedure in cartesian coordinate 
system is as below. 


* h/ *yvou 

= 0 

* u/, %oii 

= o 

dv/ *yvou 

= <2*v) 


vo.il 


/Ay 


. .. <3. 9) 
. (3. 10) 
(3. 11) 


Substitution of Eqs. 2.6.— 2.12 in Eqs. 3.9—3.11 results in 
the following equations which are valid at a solid side wall 
boundary. 


Cosor#0h/dC + Sh/&t) =0 (3.12) 

Cosa*«9u/<«; + tfu/«*T7 =0 (3.13) 

Cosot*«*v/«»C + 0 v/ 0 y i = 2 *v / At) (3.14) 


The finite difference approximation of eqs. 3-12-3.14 for 
the solid side wall shown in Fig 3.3. gives the following 
explicit equations for determining h, u and v at any node 
( i , 1 ) . 

h ( i , 1 ) = h ( i , 2) + Ar*/AC*Ch(i + 1,2> - h ( i ,2) >*Cos* ....(3.15) 
u(i,1) = u(i,2) + A-q/AC -*Ku( i + 1 ,2) - u(i,2)>*Cos« ....(3.16) 
v ( i , 1 ) = -v ( i ,2) + At)/AC»Cv( i + 1,2) - v ( i ,2)>*Cos« ...(3.15) 
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h, u and v are determined at the interior node (i, H) using 
the MacCormack scheme. These are then used in Eqs. 3.15-3.17 
to determine h, u and v at the boundary nodes (i,1). Similar 
procedure is used at all other solid side wall boundaries. 
Reflection procedures are relatively easier to implement. 
However, they have been criticised because of the introduction 
of additional conditions at the boundaries (such as Eqs. 3.9 
and 3.10) which represent symmetric conditions there by over 
specifying the problem (Moretti 1969). Roache (1972) examined 
this problem and showed that if the grid is set as shown in 
Fig 3.3. (i.e. setting the wall at y= 0.5#Ay the errors 
introduced by the reflection procedure is minimum. 

3. 3 Stability condition 

The mathematical foundations for the question of stability of 
numerical schemes are well developed only for linear systems. 
The results from linear theory are used as guidelines to 
nonlinear problems,, the justifications depending on numerical 
experiments. Governing equations for our case are nonlinear. 

MacCormack scheme is stable if the courant-FriedrichsrLewy 
(C F L ) condition is satisfied at every grid point. This 
condition for two dimensional flows is heur i s t i cal ly expressed 
by the following equations (Roache 1972) 
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Cn = -C(V -M gh >AtJ (A£ 2 +A» a ) } / (A< /At) )S 1 (3.18) 

In terms of At 

At = Cn/< v+fgF>*( 1/J ( 1/(A£ 2 ) + 1 /(Ay) 2 ) ) ) (3.19) 

In which V is resultant velocity at the grid point and Cn 
should be less than one. Above condition is based on a 
linearized form of the governing equations for one dimensional 
flows. A certain amount of numerical experimentation is 
required before choosing the actual upper limit for cn. 


3. 4 Artificial Viscosity 

MacCormack scheme is second order accurate in space and 
time and results in dispersive errors (Anderson et al 1984). 
These dispersive errors cause high-frequency oscillations near 
steep gradients. In the model presented herein, the procedure 
developed by Jamison et al. (1981) is used to dampen these 
high frequency oscillations. This procedure smoothens regions 
of large gradients while having smooth areas relatively 
undisturbed. The values of the variable at the new time 
computed by MacComack method are modified using the following 
algor ithm. 


| h(i + 1 , j )-2*h(i, j >+h( i— 1 , j ) | 

| h ( i + 1 , j ) | +| 2#h ( i , j >| +| h(i-1 , j ) j 


(3.20) 



41 


v> « h( i , j + 1 >-2*h( i , j >+h< i , j-1 > j 


j h< i , j>- I ) | « j 2 h h < i , j > j +{ h( i, j-1 > } 


At points where h(i,j— 1) does not exist 


v = |h(i,j+1)-h(i,j)| 

v y . 

j h ( i , j f - 1 ) | »- j h ( i , j ) | 

At points where h ( i , j + 1 ) does not exist 

v = j hit i , j -1 ) -h < i , j ) | 

j h ( i „ j - 1 ) | + j h ( i , j ) { 

At points where h(i— 1,j> does not exist 

j h < i + 1 , j ) ~h ( i , j ) { 

( h ( i + 1 , j ) | + 1 h ( i , j ) j 

At points where h ( i+1 , j ) does not exist 


. ( 3.21 ) 


.....'( 3 . 22 ) 


( 3.23 ) 


( 3.24 ) 


- | -h ( i , j ) +h ( i - 1 , j ) | 

| h ( i , j > j + 1 h ( i - 1 , j ) j 


tV,! „ .. - l< * max(v K . . . . ) 


( 3 . 25 ) 




....( 3 . 26 ) 


V-' . , = K rna:<(L> . . ,v . ) 

17 <l-c». n,j> . . 

Where K is a dissipation constant. The corrected values of 
•he variable / at the new time step are given by the following 


...( 3 . 27 ) 


■ q u a t i a n - 



/ < i > j > ~ /Hr j) + A + B 
Where A and B are 


(3.S8) 


4E 


A -C/(i + 1 r j)-/(i,j)J#™ _ •[/ (i , j )-/ ( i-1 , j )>^ 

B,j> ' J B,j> 

= {/<i + 1,j>-/<i,j)>*v C/(i,j>-/<i-1, j)>^ . . 

17<V,J 4o . !5) ' *•' ^7)<L, j-O . B> 

Value of variable / is the one obtained after the corrector 
step. Dissipative constant, K, is used to regulate the amount 
of dissipation. The above procedure is equivalent to adding 
second order dissipative terms to the original governing 
equations. The actual numerical eddy viscosity coefficient is 
of the order of Kv< * (A£ 2 /A t 2 ). This indicated that the 
influence of K on the results depends upon gradients in the 
flow depth as well as grid size. Equations show that influence 
of this procedure of artificial viscosity is minimal in 
smooth regions. Since yt tends to zero in such a case.Kroll and 
Jain (1987) showed that the influence of dissipation 
coefficient is reduced considerably , if the grid is refined 
and consequently a slightly larger value of K is needed to 
smooth out the oscillations. It should be noted here that the 
presence of friction term which is dissipative in nature 
tends to inhibit the oscillations. Therefore, a small value of 
K is sufficient for smoothing in such situations . 


3. S Closure 

In this chapter, MacCormack scheme for solving the 
governing equation has been presented. The boundary conditions, 
stability criteria and the artificial viscosity procedure have 
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been discussed. The principles presented here were used to 
develop a computer program for analyzing the supercritical 
flow in channel junctions. Results of this study are presented 
in the following chapter. 



CHAPTER 4 


RESULTS AND DISCUSSION 

To verify the mathematical model presented in chapter H and 3, 
computed results are compared with laboratory test data (Hager 
1989). The numerical results are also compared with the 
analytical solution for maximum flow depth (sect-1. 3) to 
determine range of parameters for which the analytical model 
correctly simulates the two-dimensional flow pattern. 



BRANCH - 2 


plan view of e xperimental arrangement. 


-Fig. 4-1 
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The experiments (Hager 1989) were performed in a 

rectangular, horizontal channel 50 cm wide. Two elements were 

inserted in the channel as shown in the Fig. 4.1, making the 

widths of the three branches equal to 99 mm each. Two 

junction angles, <5 = 22.5° and 45° were considered. The flow 

depths were controlled by two vertical, movable gates which 

were placed immediately upstream of the junction. The length 

of the downstream branch was 80 cm and after this it expanded 

again to the full channel width of 50 cm. The main test 

section was situated between the gates and the end of the 

downstream branch. Two types of test were conducted. Flow 

depths h and h,, the velocities v and v,, the angle of the 
■ u 1 u 1 

main wave front Q and the maximum flow depths h were 

max 

measured in majority of the runs. The Froude numbers covered 
the ranges 2.8< F < 16 for 6— 22.5° and 3.3 < F < 8.3 for 6= 
45°. The entire three-dimensional velocity field and the flow 


surfaces 

were recorded for three 

selected 

runs . 

These 

selected 

runs were made for <5 = 22.5° 

. A point 

gauge 

<± 0.1 

mm ) an d a 

propeller velocity meter (± 

3 cm/s) 

were 

used to 


take measurements. 


4. i Comparison For 2-D Watar Surf aca Prof ila 

The inflow conditions for which the entire water surface 
profile was measured are given in Table 4.1 






Fig. 4 


J * 


(b ) Measured 


.5 Water level Contour for Run 1 


oot 
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Tabl« 4.1 

Experimental details of the three selected runs 


Run # 

F 

u 

F 1 

h u (cm) 

h T (cm) 

I 

5.25 

6.2 

2.65 

1.98 

2 

4.5 

4.5 

2.65 

2.65 

3 

j 

6.2 

5.25 

2.65 

3.58 


The mathematical model was run for the above test 

conditions using a numerical grid of AC = A 6 - 3.5 cm. This 

corresponds to 9 nodes in the main channel. Choice of grid 

size was based on balancing computational speed with the 
desired accuracy. The computational time step. At was 
controlled by the courant number and grid size. As mentionned 
earlier the upper limit for the Courant number to achieve 
numerical stability required some amount of numerical 
experimentation. The Courant number for the runs ranged from 
0.65-0.85. The dissipation coeff. K, was determined by trial 
and error. As discussed in chapter 3, the influence of K upon 
the results depends on gradients in flow depths as well as on 
grid size. K should be as small as possible and at the same 
time it should reducethe higher order numerical oscillations. 
The K value for the above runs was 0.9. The roughness coeff. 
is assumed to be zero since no information on this was 

CENTUM- 1 

— ,1IZ5S£ 

isc.No. 




TABLE 4. 2 ; <5- 22. S & 


# 

h 

u 

h l 

F 

u 

F 1 

h 

a 

h n 

h e 

1 

2.65 

2.65 

4.26 

4.18 

8.84 

9.5 

14.3 

2 

2 * 68 

2.68 

4.61 

4.47 

9.66 

10.2 

16.0 

3 

2.67 

2.67 

5.29 

5.30 

11.3 

11.6 

19.8 

4 

2.65 

1.76 

3.94 

4.72 

7.58 

7.89 

9.9 

5 

2.65 

1 .77 

4.37 

5.23 

8.41 

7.88 

11.2 

6 

2.65 

1 .78 

4.80 

5.72 

9.23 

8.36 

11.9 

7 

2.65 

1 .78 

5 « 22 

6.22 

10.0 

9.33 

13.6 

8 

2.71 

1 .33 

3.78 

5.48 

7.13 

7.0 

8.5 

9 

2.72 

1.33 

4 » 65 

6.53 

8 » 62 

7.5 

9.9 

10 

2.72 

1.34 

5.15 

7.17 

9.50 

8.25 

11.3 

11 

2.74 

0.67 

5.07 

10.1 

8.07 

7.0 

6.7 

12 

1.85 

2.55 

3.78 

3.12 

5.71 

7.43 

9.2 

13 

1.85 

2.60 

4.15 

3.60 

6.45 

8.30 

11.1 

14 

1 .87 

2.53 

5.09 

4.37 

8.0 

9.44 

13.2 

15 

1.32 

2.68 

6.78 

4.70 

8.34 

10.2 

14.7 

16 

0.65 

2.68 

9.22 

4.52 

6.68 

9.79 

15.4 

17 

0.35 

2.54 

14.8 

5.42 

7.2 

11.2 

18.4 

18 

0.35 

0.34 

13.1 

15.1 

4.86 

6 . 80 

8.1 

19 

0.75 

0.37 

6.05 

8.78 

3.13 

3.2 

3.8 

20 

0.75 

0.37 

9.14 

13.1 

4.79 

6.8 

6.2 

21 

0.73 

0.73 

4.76 

4.76 

2.72 

2.85 

3.7 

22 

0.73 

0.73 

7.18 

7. 18 

4.27 

4.77 

5.6 

23 

0.72 

0.72 

9.67 

9.67 

5.9 

7.76 

10.2 

24 

0.32 

1.31 

10.5 

5. 1 

3.87 

5.39 

6.6 


h = maxm. flow depth 

Subscripts a, n, e denotes analytical, numerical 
and experimental respectively. 
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TABLE 4.3 } 6m 45. 


# 

h u 

h l 

F 

u 

F 1 

“h 

a 

h n 

h e 

-1 

1.24 

1 .28 

4.64 

4.37 

7.33 

9.6 

9.2 

2 

1.25 

1 .29 

5.80 

5.45 

9.96 

13.6 

12.3 

3 

1.28 

1 .32 

6.72 

6.05 

12.2 

17.2 

15.6 

4 

2.48 

2.50 

4.08 

3.92 

12.3 

16.9 

15.7 

5 

2.46 

2.52 

4.48 

4.32 

13.9 

18.7 

17.4 

6 

1.83 

2.46 

4.96 

4.17 

12.4 

17.3 

15. 1 

7 

1.79 

2.53 

5.54 

4.58 

14.2 

19.9 

17.2 

8 

1.22 

2.46 

6.30 

4.52 

12.3 

17.4 

14.7 

9 

1.23 

2.47 

7.00 

4.82 

14.2 

19.9 

18.4 

10 

2.43 

1.85 

3.79 

4* m 34* 

10.5 

13.9 

12.9 

11 

2.46 

1.85 

4.28 

4.81 

12.5 

16.4 

15.6 

12 

2.47 

1.84 

4.71 

5.25 

14.2 

18.4 

18.4 

13 

2.51 

1.30 

3.79 

5.10 

10.2 

12.7 

12.9 

14 

2.48 

1.28 

4.42 

5.93 

12.2 

15.7 

15.0 

15 

2.47 

1.30 

4.94 

6.55 

14.1 

18.6 

18.0 

16 

2.45 

0.68 

3.30 

6.58 

7.7 

9.64 

08.6 

17 

2.48 

0.70 

3.89 

6 . 98 

9.2 

11.6 

10.0 

18 

0.70 

0.75 

5.53 

5. 16 

5.30 

6.90 

7.7 

19 

0.71 

0.76 

8.11 

7.40 

8.89 

15.0 

13.0 

20 

0.46 

0.50 

5.41 

4.74 

3.32 

4.27 

4.5 

21 

0.47 

0.51 

6.70 

6.12 

4.59 

6 . 64 

6.2 

22 

0.49 

0.53 

7.07 

6 . 53 

5.16 

7.27 

7.6 

23 

0.51 

0.56 

7. 15 

6.53 

5.46 

7.25 

8.5 


h « maxm. flow depth 

Subscripts a, n, e denotes analytical, numerical 
and experimental respectively. 
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